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Abstract:  We  present  a  new  hybrid  explicit/implicit  solvent  method  for  dynamics  simulations  of  macromolecular 
systems.  The  method  models  explicitly  the  hydration  of  the  solute  by  either  a  layer  or  sphere  of  water  molecules,  and 
the  generalized  Born  (GB)  theory  is  used  to  treat  the  bulk  continuum  solvent  outside  the  explicit  simulation  volume.  To 
reduce  the  computational  cost,  we  implemented  a  multigrid  method  for  evaluating  the  pairwise  electrostatic  and  GB 
terms.  It  is  shown  that  for  typical  ion  and  protein  simulations  our  method  achieves  similar  equilibrium  and  dynamical 
observables  as  the  conventional  particle  mesh  Ewald  (PME)  method.  Simulation  timings  are  reported,  which  indicate 
that  the  hybrid  method  is  much  faster  than  PME,  primarily  due  to  a  significant  reduction  in  the  number  of  explicit  water 
molecules  required  to  model  hydration  effects. 
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Introduction 

In  the  computer  simulation  of  biomolecules,  one  often  has  the 
choice  of  treating  the  aqueous  solvent  as  a  collection  of  explicit 
water  molecules1,2  or  as  an  implicit  solvent  model  such  as  the 
generalized  Born  (GB)  theory.3  Implicit  solvent  simulations  are 
typically  faster  than  fully  microscopic  evaluations  and  easier  to 
interpret  as  the  water  degrees  of  freedom  are  absent.  Because  of 
their  mean-field  approximation,  implicit  solvent  models  are  of 
lower  resolution  and  have  been  known  to  blur  the  potential  energy 
landscape  of  a  protein,4  and  cause  structural  distortions.5  More¬ 
over,  implicit  models  lack  local  hydrogen  bond  interactions  be¬ 
tween  the  solute  and  solvent  that  can  lead  to  incorrect  preferences 
for  secondary  structural  motifs.6 

Alternatively,  explicit  solvent  methods  offer  a  more  detailed 
and  accurate  description  of  a  macromolecular  system.  Yet,  these 
methods  often  expend  a  majority  of  their  computational  effort  on 
the  simulation  of  water  molecules  rather  than  the  solute.  Further¬ 
more,  certain  periodic  boundary  methods,  such  as  particle  mesh 
Ewald  (PME),  evoke  artificial  real- virtual  solute  interactions  that 
can  be  problematic  when  determining  thermodynamic  quantities 
such  as  free  energies.7-9 

Methods  that  combine  explicit  and  implicit  solvent  ap¬ 
proaches2,10-12  seek  to  gain  the  salient  features  of  both  ap¬ 


proaches;  namely,  accuracy  and  speed.  There  are  two  classes  of 
hybrid  explicit/implicit  solvent  schemes.  In  one  class,  known  as 
reaction  field  methods,12-14  the  electrostatic  interactions  between 
molecules  are  radially  truncated,  and  a  dielectric  continuum  be¬ 
yond  the  cutoff  is  employed  to  estimate  the  dielectric  screening  of 
the  neglected  interactions.  An  advantage  of  this  method  is  that 
periodic  boundary  conditions  can  be  employed,  thus  allowing 
unrestricted  movement  of  water  molecules.  A  further  benefit  is  the 
use  of  rather  short  nonbonded  cutoffs  (—10  A).  A  disadvantage  is 
that  the  procedure  assumes  that  a  dielectric  of  80  always  lies 
outside  the  cutoff  region,  which  is  an  inappropriate  assumption 
when  modeling  large  solutes  such  as  biomolecules.  An  additional 
problem  is  that  the  reaction  field  method  uses  the  same  type  of 
periodic  box  as  a  conventional  Ewald  method  and  thus  requires  the 
same  number  of  water  molecules  to  solvate  a  system. 

A  second  type  of  hybrid  solvent  approach,  which  we  will 
designate  as  cluster  methods,2,10,15,16  utilizes  a  simulation  volume, 
often  spherical  in  shape,  surrounded  by  a  dielectric  continuum. 
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Unlike  procedures  using  periodic  boundary  conditions  (PBC), 
cluster  methods  do  not  allow  free  movement  of  water  molecules 
near  the  simulation  surface  and  thus  suffer  from  surface  artifacts17 
that  can  be  partially  corrected.2,10,15,18  Cluster  methods  are  not 
plagued  by  the  artificial  interactions  between  real  and  virtual  solute 
atoms  found  in  periodic  boundary  Ewald  methods.7  Furthermore, 
cluster  methods  have  been  shown  to  provide  convergent  observ¬ 
ables  with  increasing  simulation  radius.7,10  Unlike  reaction  field 
methods,  cluster  approaches  cannot  be  reliably  used  with  short 
nonbonded  cutoffs;  however,  long-range  electrostatic  approxima¬ 
tions  that  greatly  reduce  computational  effort  have  been  devised.8 
The  primary  savings  in  finite  cluster  simulations  over  conventional 
periodic  boundary  Ewald  approaches  is  achieved  because  fewer 
numbers  of  water  molecules  are  required.2,10 

Still,  there  are  many  difficulties  with  cluster  approaches.  First, 
a  spherical  boundary  is  often  used  because  there  is  an  analytical 
solution  to  the  reaction  field  of  a  spherical  dielectric.10,19  A  spher¬ 
ical  shell  of  water  molecules  can  be,  however,  quite  wasteful  for 
highly  nonspherical  solutes.10,20  Simulation  volumes  that  do  not 
have  simple  geometric  shapes  require  numerical  solution  of  the 
Poisson  equation  (or  Langevin  equation),  which  is  a  comparatively 
slower  process.20  A  second  problem  with  cluster  approaches  when 
used  with  an  analytical  Poisson  formula,  is  that  the  reaction  field 
potential  approaches  infinity  at  the  boundary.  This  is  due  to  the 
incorrect  assumption  that  individual  atoms  do  not  contribute  their 
own  dielectric  volume  as  they  pierce  through  the  spherical  bound¬ 
ary. 

There  are  many  cluster  approaches  that  employ  irregular  sim¬ 
ulation  volumes  that  follow  the  shape  of  the  biological  molecule  to 
reduce  the  number  of  explicit  solvent  molecules.  For  example, 
Beglov  and  Roux21  introduced  a  method  that  uses  a  sum-over- 
spheres  simulation  volume  and  maintains  constant  pressure22,23 
through  dynamic  boundaries.  Unfortunately,  this  approach  and 
most  other  irregularly  shaped  cluster  approaches  lack  a  reaction- 
field  treatment.  The  one  exception  is  the  SAPHYR  method  of 
Lounnas  et  al.,11  which  uses  an  image  dipole  formulation  to  treat 
the  reaction  field  of  each  water  molecule.  Nonetheless,  it  does  not 
appear  that  the  SAPHYR  method  incorporates  a  reaction  field  for 
the  solute  atoms.  Perhaps  this  is  because  the  position  of  the  image 
dipole  is  undefined  for  a  solute  atom  as  the  atom  sits  at  the  center 
of  the  sphere  it  carves  out. 

In  this  work,  we  present  a  method  that  also  uses  a  generalized 
sum-over- spheres  boundary  but  incorporates  a  continuum  reaction 
field  via  GB  theory.  The  GB  solvent  model  has  been  shown  to  be 
a  reasonable  alternative  to  the  solution  of  the  Poisson  equation, 
obtaining  an  average  absolute  error  of  —1%  for  a  large  set  of 
proteins.24,25  GB  solvation  energies  and  forces,  unlike  Poisson 
theory,  only  require  a  single  noniterative  solution,  and  hence  are 
much  faster  than  Poisson  methods  for  arbitrary  dielectric  volumes. 
In  our  approach,  the  dielectric  boundary  is  held  fixed,  which 
further  simplifies  the  problem.  With  a  fixed  dielectric  boundary, 
the  Born  radius  becomes  a  scalar  function  that  can  be  stored  on  a 
grid.  This  approximation  enables  atomic  Born  radii  to  be  obtained 
expediently  by  interpolation. 

With  the  time-consuming  Born  radii  integration  step  removed, 
the  bottlenecks  in  our  method  become  the  pairwise  electrostatic 
and  GB  terms.  Use  of  a  simple  nonbonded  cutoff  in  explicit  water 
simulations  is  well  known  to  degrade  accuracy.7  Approximate 


treatments  of  the  long-range  electrostatic  component  include  the 
extended  electrostatics  approach,26  the  local  reaction  field 
method,8  and  the  fast  multipole  method.27  Nonetheless,  it  is  not 
obvious  how  to  modify  these  methods  to  handle  GB  terms  because 
they  are  all  based  on  mathematical  manipulations  of  the  Coulom- 
bic  kernel.  Fortunately,  the  multigrid  approach28,29  is  an  alterna¬ 
tive  strategy  for  reducing  computational  cost  for  pairwise  interac¬ 
tions  that  does  not  require  a  Coulombic  kernel.  The  multigrid 
approximation  to  pairwise  point  charge  interactions  has  linear 
scaling  computational  behavior  with  respect  to  system  size  and 
early  crossover  with  respect  to  exact  calculations.29  Furthermore, 
the  accuracy  of  the  multigrid  approach  is  tunable  and  within  the 
error  tolerances  required  for  molecular  dynamics  simulations.29  In 
this  work,  we  modified  the  multigrid  approach  of  Skeel  et  al.29  to 
approximate  the  long-range  portion  of  the  pairwise  GB  term. 

In  the  remainder  of  this  work,  we  present  the  details  of  our 
hybrid  solvent  method  and  multigrid  algorithm.  Then,  we  present 
results  showing  the  computational  benefits  of  the  multigrid  proce¬ 
dure  over  a  standard  cutoff  approach.  Next,  we  compare  dynamical 
and  structural  results15,30  of  several  model  simulations  between 
our  procedure  and  the  conventional  PME  explicit  solvent  method. 
Finally,  we  envisage  the  types  of  applications  where  the  hybrid 
method  may  have  an  advantage  over  conventional  PBC  ap¬ 
proaches. 


Theory  and  Methods 

Definition  of  Hybrid  Solvent  Potential 

The  goal  of  a  hybrid  solvent  approach  is  to  partition  solvent 
degrees  of  freedom  into  explicit  and  implicit  regions.10,31,32  In  our 
case,  a  simulation  volume  is  defined  as  a  region  in  space  that 
contains  the  solute  and  explicit  water  molecules.  Outside  the 
simulation  volume,  an  implicit  solvent  is  described  by  a  high 
dielectric  continuum.  The  dielectric  continuum  can  be  thought  of 
as  thermal  reservoir  of  infinitesimal  dipoles  that  interact  with  the 
solute,  explicit  water  atomic  charges,  and  themselves.  Figure  1 
schematically  illustrates  the  system  setup  for  a  hybrid-layer  sim¬ 
ulation. 

Previous  reaction  field  approaches  have  often  been  limited  to 
spherical  volumes,  as  one  of  the  few  analytic  solutions  to  the 
Poisson  equation  is  for  a  spherical  boundary.10  GB  theory  offers  a 
fast,  yet  approximate,  analytical  treatment  of  the  dielectric  contin¬ 
uum  that  is  not  limited  to  a  spherical  boundary.  However,  we  have 
chosen  a  rather  simple  dielectric  volume  that  is  composed  of  a 
transformed  sum-over-spheres.  Combining  spheres  of  radius  w 
around  each  heavy  atom  is  nearly  identical  to  describing  a  layer  of 
width,  w,  surrounding  the  entire  solute.  To  assure  continuity  at  the 
intersection  of  spheres  centered  on  different  atoms,  we  use  a 
volume  function  similar  to  the  one  proposed  by  Im  et  al.33  If  we 
define  a  volume  function,  V(r),  such  that  V  =  1  is  defined  as  inside 
and  V  =  0  is  defined  as  outside,  the  interlocking-sphere  volume 
has  the  following  form: 


V(r) 


i-n 


i  - 


(i) 
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Figure  1.  Schematic  of  the  hybrid  explicit/implicit  solvent  approach. 
The  simulation  volume  consists  of  a  solute  immersed  in  a  layer  of 
water  molecules.  The  gray  outer  region  corresponds  to  the  bulk  di¬ 
electric  continuum. 


where 


f  1  s  <  0 

t»(s)  =  |  1  —  6s5  +  15s4  —  10s3  0  <  s  <  1,  (2) 

[  0  s  >  1 

xt  is  the  position  of  the  center  i,  and  w°  and  w1  are  the  starting  and 
stopping  points,  respectively,  of  the  atomic  function  tail.  Given 
this  volume  function,  it  is  now  possible  to  define  the  reaction  field 
and  repulsive  boundary  potential;  however,  to  increase  flexibility 
in  parameterization,  the  values  of  w°  and  w1  are  different  for  the 
two  energy  terms  (see  below). 

The  reaction  field  is  defined  using  GB  theory.  Because  the 
simulation  and  dielectric  volumes  are  assumed  fixed,  Born  radii 
are  precalculated  on  a  three-dimensional  grid  with  a  cell  size  of  2 
A.  During  the  simulation,  the  Born  radius  for  each  atom  is  ob¬ 
tained  by  cubic  interpolation  (see  below)  of  the  Born  radii  grid. 
The  Born  radii,  am,  of  each  grid  point,  m  is  defined  in  terms  of  the 
volume  function,  V(r ): 


where  R  is  an  asymptotic  radius  set  to  1.5  A,  and  xm  is  the  position 
of  grid  point,  m.  After  am  is  calculated  via  eq.  (3)  it  is  linearly 
transformed,  a'm  =  1.0649am  —  0.316.  This  transformation  was 
derived  empirically  by  fitting  the  GB  self  energies  of  a  1000 
random  points  inside  of  a  sphere  to  the  analytical  Poisson  results.10 
Equation  (3)  is  evaluated  via  a  numerical  integration  algorithm  that 


uses  a  26-point  Lebedev  angular  grid  and  17  radial  points.25 
Furthermore,  eq.  (3)  has  been  shown  elsewhere25,34  to  provide  an 
accurate  estimate  of  Born  radii  compared  to  the  Poisson  method, 
with  correlation  coefficients  of  0.99. 

Given  Born  radii,  the  total  GB  energy  is  calculated  via  a 
formula  similar  to  the  one  originally  proposed  by  Still  et  al.,3 


^GB 


_ m _ 

i 

rij  +  2  («/  +  «/)exp[-2 rijlioti  +  a,-)] 


(4) 


where  qL  are  atomic  charges  and  rtj  is  distance  between  atoms  i  and 
j,  esolv  is  the  dielectric  constant  of  the  solvent  (for  this  work, 
£soiv  =  80)  and  k  =  -166.0.  We  used  this  alternative  formula 
because  it  appeared  to  provide  better  agreement  with  the  exact 
results  of  a  spherical  system  (results  not  shown).  Also,  eq.  (4) 
reduces  the  number  of  square  root  evaluations  by  a  factor  of  2  vs. 
the  original  formula.  Nonetheless,  a  variety  of  alternative  GB 
formulas  are  possible  that  provide  similar  results  as  the  original 
Still  formula.35,36  The  wall  potential  in  our  formulation  is  also 
defined  in  terms  of  the  volume  function  [eq.  (1)]: 


-^wallC^*) 


1 

£S£f 1  +  V(r) 


1 


/-max' 

^wall 


+  1  ’ 


(5) 


such  that  V=  1  maps  to  Ewall(r)  =  0  and  V  =  0  maps  to  £wall(r)  « 
iiwan-  Because  we  use  a  finite  value  for  F™aP,  the  wall  potential  is 
actually  penetrable.  This  turns  out  to  be  useful  because  a  system 
that  is  started  with  too  many  water  molecules  will  actually  spill  the 
excess  water  molecules  beyond  the  simulation  volume,  thus  re¬ 
lieving  some  pressure  on  the  system.  We  did  not  incorporate  an 
attractive  van  der  Waals  (vdW)  potential  because  initial  tests 
indicated  that  such  a  potential  had  a  relatively  small  impact  on 
radial  distribution  functions.  A  more  important  factor  was  the 
proper  placement  of  the  dielectric  boundary  for  the  reaction  field. 
We  also  did  not  incorporate  energy  terms  that  would  improve  the 
angular  isotropy  of  water  molecules  near  the  cluster  surface.2,18 
Such  correction  terms  are  supposed  to  reduce  artificial  forces  that 
arise  from  surface  dipoles.  However,  one  key  study  indicated  that 
these  terms  fall  short  in  reducing  surface  dipole  artifacts  in  the 
application  of  calculating  ionic  charging  free  energies  in  solu¬ 
tion.17  Furthermore,  implementing  angular  correction  terms  for 
nonspherical  boundaries  is  nontrivial  and  beyond  the  scope  of  this 
work. 

Although  we  present  the  hybrid  method  as  a  means  to  incor¬ 
porate  layer  of  explicit  water  molecules  around  a  solute,  it  can  also 
be  used  to  simulate  a  spherical  cluster  of  water  molecules.  First,  if 
a  solute  is  present,  its  center  of  mass  is  translated  to  the  origin. 
Next,  a  fixed  dummy  atom  is  defined  at  the  origin.  This  dummy 
atom  is  the  only  atom  used  to  define  the  simulation  volume.  The 
simulation  sphere  radius,  ^sphere,  is  calculated  as  a  sum  of  the 
specified  layer  width  w  and  the  maximal  distance  between  a  solute 
atom  and  the  origin.  This  simple  protocol  could  be  improved, 
because  it  does  not  provide  the  optimally  smallest  simulation 
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Figure  2.  The  interaction  kernels  used  in  the  multigrid  implementation  of 
the  Coulomb  electrostatic  term.  In  this  example,  the  cutoff  value,  a,  is  8  A. 
Dashed  line:  Coulomb  kernel,  1/r,  solid  line:  local  kernel,  Klocal(r),  dot- 
dashed  line:  soft  kernel,  Ksoft(r). 


sphere  that  can  encompass  the  solute  and  still  provide  at  minimum 
a  width,  w,  of  water  molecules.37 


and  combine  with  the  local  interactions  in  step  1  to  obtain  the 
energy  and  forces  (L  =  1  — >  L  =  0). 


Multigrid  of  Pairwise  Interactions 

Having  simplified  the  normally  expensive  Born  radii  calculation, 
the  Coulomb  and  GB  pairwise  interactions  become  the  most  com¬ 
putationally  expensive  terms  in  the  hybrid  solvent  method.  Thus, 
we  implemented  a  multigrid  approach  to  approximate  the  long- 
range  component  of  the  electrostatic  and  GB  interactions,  and 
hence  greatly  reduce  the  computational  cost.  The  multigrid  method 
for  Coulomb  electrostatics  has  been  developed  by  Skeel  et  al.29 
and  Sandak.28  We  adapted  the  multigrid  method  of  Skeel  et  al.29 
to  incorporate  the  GB  pairwise  screening  formula  [eq.  (4)]. 

The  basic  principle  of  a  multigrid  electrostatics  procedure  is  to 
cast  pairwise  charge  interactions  onto  grids  through  interpolation 
and  treat  various  interaction  distances  with  different  degrees  of 
approximation  or  grid  resolutions  (indexed  as  L).  The  multigrid 
algorithm  is  sketched  as  follows: 


To  separate  the  ranges  in  step  4  precisely,  the  interaction  kernel 
for  each  grid  level  L,  KL(r ),  is  split  into  two  components:  a  soft 
function  and  a  local  function, 

K\r )  =  KLsoft(r)  +  ^focal(r).  (6) 

The  soft  function  for  the  L  =  0  interactions  must  slowly  vary 
because  it  is  the  basis  of  interaction  for  the  L  =  1  grid.  For  a 
Coulombic  potential,  K°(r)  =  1/r,  the  soft  function  is  chosen  such 
that  £?ocal(r)  goes  to  zero  beyond  cutoff  a.  In  the  original  article  of 
Skeel  et  al.,29  several  soft  functions  are  presented  and  tested. 
Nonetheless,  we  desired  to  reduce  the  variance  of  the  soft  function 
completely  in  the  domain,  r  <  all ,  so  that  self-interaction  and 
exclusion  interaction  artifacts  could  be  reduced  or  eliminated. 
Therefore,  we  created  a  new  piecewise  soft  potential  that  is  con¬ 
tinuous  up  to  second  derivatives  at  r  =  all  and  r  =  a; 


1.  Perform  explicit  local  interactions  of  atomic  charges,  when  r  < 
a,  where  a  is  user-defined  local  cutoff  (L  =  0). 

2.  Interpolate  atomic  charges  onto  the  finest  resolution  charge  grid 
(L  =  0  -»  L  =  1). 

3.  Create  a  hierarchy  of  grids  by  interpolating  each  charge  grid  to 
the  charge  grid  with  the  next  coarsest  resolution  (L  =  n  — >  L  = 
n  +  1). 

4.  For  each  individual  charge  grid,  L,  build  a  potential  grid  such 
that  each  grid  point  sees  only  charges  less  than  a  distance  lLa 
away. 

5.  Interpolate  potential  grids  from  coarse  to  fine  resolution  (L  = 
n  +  1  — »  L  =  n). 

6.  Contract  atomic  charges  with  the  finest  resolution  potential  grid 


31 

14a  ’ 


a\3  31  a 

2/  +  24a  ’  2<r<0 


(7) 


The  local  function  for  L  =  0,  by  definition,  is  the  difference  of 
the  original  kernel  and  the  soft  function, 


^Hocal(r) 


(8) 
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The  interaction  kernels  for  L  =  0  are  depicted  in  Figure  2. 
Analogously,  for  each  grid  level,  L,  except  the  coarsest,  L  =  Lmax, 
the  interaction  kernel  is  composed  of  a  soft  function  where  a  in  eq. 
(7)  is  replaced  by  lL~l  a ,  minus  a  softer  function,  where  a  is 
substituted  with  2 La.  For  the  coarsest  grid,  Lmax,  the  interaction 
kernel  is  simply  the  softest  function  corresponding  to  a  in  eq.  (7) 
being  replaced  by  2Lmax~l  a.  In  this  work,  Lmax  was  set  to  3.  Larger 
values  of  Lmax  might  be  required  for  system  sizes  larger  than  the 
ones  treated  here. 

Multigrid  steps  2,  3,  5,  and  6  involving  cubic  interpolation 
utilize  six-dimensional  basis  functions  of  the  form:29 


is  also  split  into  local  and  soft  terms,  /local  and  Jsoft: 

JL(r ip  oii,  aj )  =  JjocXrip  aj)  +  Jift(ri},  a„  a,-).  (11) 

Some  experimentation  determined  that  a  reasonable  splitting  in¬ 
volves  setting  rtj  equal  to  the  corresponding  inverse  of  the  respec¬ 
tive  Coulomb  terms,  K^ocal  and  K^oft: 

jLa\(rip  Oth  aj )  =  7(K,cai]"\  <*;,  aj) 

Jiffrij,  OLi,  aj)  =  /(Koft]-‘,  a„  aj).  (12) 


ya,  Za,  xb,  yb,  Zb)  =  <l> 


(9) 


where 


Furthermore,  Born  radii  are  evaluated  at  every  multigrid  point. 
Born  radii  at  the  lowest  level  multigrid  (L  =  1)  are  obtained  via 
bilinear  interpolation  of  the  same  Born  radii  grid  used  to  determine 
atomic  Born  radii.  Coarser  grids  (L  >  1)  are  obtained  by  cubic 
interpolation  of  their  next  finest  resolution  Born  radii  grid  (L  —  1). 


<*>(/>) 


(i-LH)(i  + 


>1 


D(2  -  H)2 


Ip  I  s  1 

1  <  [p|  <  2  , 

Ip  I  —  2 


(10) 


where  (jta,  ya,  za)  and  (xb,  yb,  zh)  are  3D  coordinates  corresponding 
to  the  source,  a,  and  destination,  b ,  of  a  charge  or  potential,  L  is  the 
level  of  the  destination  grid,  and  hL  is  the  side  length  of  a 
destination  grid  cube.  For  example,  in  step  2,  (xa,  ya,  za)  corre¬ 
sponds  to  the  atomic  coordinates,  and  (xb,  yb9  zb)  refers  to  the  L  = 
grid  coordinates.  In  step  3,  (xa,  ya,  za)  signifies  the  L  =  n  grid 
coordinates  and  (xb,  yb,  zb)  corresponds  to  the  L  =  n  +  1  grid 
coordinates.  In  this  work,  the  highest  resolution  grid  was  assigned 
a  cell  size,  hL  =  4  A.  This  value  was  chosen  based  on  a  compro¬ 
mise  between  speed  and  accuracy.  A  detailed  analysis  using  var¬ 
ious  values  of  hu  several  interpolation  formulae,  and  alternative 
soft  functions  is  presented  in  Skeel  et  al.29 

The  multigrid  procedure  is  one  of  the  simplest  algorithms  for 
reducing  the  computational  cost  of  long-range  electrostatic  inter¬ 
actions;  however,  certain  artifacts  associated  with  casting  an 
atomic  problem  onto  a  grid  are  to  be  anticipated.  For  instance,  the 
multigrid  algorithm  erroneously  includes  the  excluded  self,  1-2, 
1-3,  and  1-4  Coulomb  interactions.  These  interactions  are  approx¬ 
imately  subtracted  out  from  the  energy  by  simply  including  the 
explicit  short-range  terms  without  the  Coulomb  kernel  term  [eq. 
(8)].  Nevertheless,  grid  artifacts  remain  unless  other  strategies  are 
employed.28,29  For  example,  in  a  test  calculation  that  we  con¬ 
ducted,  the  Coulombic  self-energy  of  an  isolated  sodium  ion 
ranged  from  1.5  to  3.8  kcal/mol,  depending  on  its  position  relative 
to  the  grid.  We  did  not  choose  to  remove  grid  artifacts  in  this  work; 
nevertheless,  understanding  the  effects  of  these  artifacts  and  re¬ 
moving  them  should  be  the  topic  of  future  work. 

To  extend  the  multigrid  procedure  to  handle  the  GB  solvent 
model,  the  kernel  from  eq.  (4), 

2=1  +  2  (<*,■  +  a,)exp[— 2 rtJ(at  +  a,)]! 


Hybrid  Simulation  Protocol 

The  hybrid  solvent  method  described  in  this  work  was  incorpo¬ 
rated  into  the  CHARMM  program  (version  c30bl).38  Molecular 
dynamics  simulations  in  this  work  were  performed  with  the 
PARAM22  empirical  forcefield  potential.39  The  equations  of  mo¬ 
tion  were  integrated  with  a  Langevin  dynamics  algorithm  at  con¬ 
stant  temperature  (298  K)  using  a  Langevin  friction  constant  of  1 
ps-1  and  an  integration  time  step  of  2  fs.  Covalent  bonds  between 
heavy  atoms  and  hydrogens  were  constrained  using  the  SHAKE 
algorithm.40  In  the  dynamics  simulations,  the  local  cutoff  for  the 
multigrid  procedure,  a ,  was  set  to  8  A.  The  vdW  terms  were 
truncated  with  a  simple  cutoff  that  switched  off  from  7  to  8  A.  Pure 
solvent  and  single  ion  simulations  (Na+,  Cl-)  were  run  for  a  total 
simulation  time  of  1  ns. 

Two  hybrid-based  molecular  dynamics  simulations  of  a  protein 
were  performed  in  this  work.  One  simulation  involved  a  10- A 
layer  (w  =  10  A)  of  water  molecules  surrounding  the  protein.  The 
other  simulation  involved  a  spherical  cluster  of  water  molecules 
with  at  least  10- A  coverage  for  every  atom  on  the  protein  surface. 
We  also  ran  a  layer  calculation  that  excluded  the  entire  GB 
component,  but  otherwise  was  identical  to  the  hybrid  layer  simu¬ 
lation.  The  choice  of  a  10- A  layer  was  based  on  the  observation 
that  charging  free  energies  of  ions  approximately  converge  at  that 
size  (results  shown  below).7,41  Production  simulations  were  run 
for  2  ns  with  structures  saved  every  picosecond.  The  protein  used 
in  this  work  was  the  62-residue  B1  domain  of  protein  L  (PDB: 
1HZ6).42  The  protein  and  water  molecules  were  first  minimized 
for  100  steps  using  the  adopted-basis  Newton-Raphson  method  to 
remove  bad  vdW  contacts.  This  was  followed  by  allowing  the 
water  molecules  to  equilibrate  for  10  ps  with  the  solute  heavy 
atoms  held  fixed. 

The  general  procedure  for  preparing  a  hybrid  solvent  calcula¬ 
tion  consists  of  two  steps:  (1)  determining  of  the  number  of  water 
molecules  needed  to  fill  the  simulation  volume,  and  (2)  carving  out 
these  water  molecules  in  the  shape  of  a  layer  or  sphere  from  a  large 
block.  The  target  number  of  water  molecules  is  calculated  as  a 
product  of  the  water  simulation  volume  (WSV)  and  a  specified 
bulk  water  density  of  0.0334  A-3.  The  WSV  is  equal  to  the 
difference  between  the  simulation  volume  and  the  standard  solute 
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Table  1.  Accuracy  of  Multigrid  (MG)  and  Standard  Cutoff  (STD)  Approaches  for  the  Model  System  of  a 
Sodium  Ion  Embedded  in  a  17-A  Sphere  of  Water  Molecules. 


Method 

Cutoff  start  (A) 

Cutoff  stop  (A) 

%  EE 

%  AFE 

%  MFE 

MG 

8 

8 

0.24 

3.4 

12 

MG 

10 

10 

0.053 

1.7 

8.7 

MG 

12 

12 

0.014 

0.97 

3.9 

STD 

8 

10 

0.13 

49 

220 

STD 

10 

12 

0.50 

35 

160 

STD 

14 

16 

3.0 

19 

94 

STD 

18 

20 

0.75 

8.4 

36 

STD 

22 

24 

0.51 

3.1 

15 

STD 

26 

28 

0.0018 

0.78 

4.4 

Definitions  of  error  measures  are  described  in  the  text.  Energies  and  forces  consist  of  the  electrostatic  and  GB  terms 
only. 


vdW  volume.  For  the  layer  simulations,  the  simulation  space  is 
obtained  by  setting  all  of  the  heavy  atomic  radii  to  w.  The  coor 
volume  command  in  CHARMM  was  used  to  numerically  calculate 
the  simulation  volume  and  the  solute  vdW  volume.  For  the  spher¬ 
ical  boundary  simulations,  the  simulation  volume  is  equal  to 
47t/3  Sphere-  With  a  target  number  of  water  molecules  determined, 
a  large  cubic  box  of  water  molecules  is  overlaid  onto  the  solute. 
First,  the  water  molecules  that  are  less  than  3.1  A  from  a  solute 
heavy  atom  are  deleted.  Next,  for  the  layer  approach,  an  iterative 
procedure  is  performed  in  which  water  molecules  beyond  a  certain 
cutoff  distance  from  all  heavy  atoms  of  the  solute  are  deleted.  The 
procedure  repeats  deletion  with  decreasing  cutoffs  until  the  num¬ 
ber  of  water  molecules  remaining  is  less  than  or  equal  to  the  target 
number  of  water  molecules.  The  spherical  simulations  utilize  in¬ 
stead  a  slowly  decreasing  radius.  Often,  in  actual  system  setups, 
some  water  molecules  will  be  initially  positioned  outside  the  wall 
potential.  Due  to  the  finite  nature  of  the  wall  potential,  these  water 
molecules  will  drift  away  during  a  simulation.  Thus,  we  use  the 
MMFP  module  in  CHARMM  to  place  a  large  spherical  boundary 
potential  around  the  simulation  volume  to  prevent  the  unbounded 
drift  of  water  molecules  that  initially  escape  from  the  smaller  finite 
boundary  of  eq.  (5). 

The  placements  of  the  dielectric  boundary  and  wall  potential 
were  empirically  parameterized  to  achieve  a  reasonable  radial 
distribution  function  for  a  spherical  droplet  of  bulk  water  with  a 
radius  of  w.  The  wall  potential  was  placed  between  w0(wall)  = 
w  -  0.6  A  and  w,  (wall)  =  w  +  0.4  A  with  a  height, 
arbitrarily  set  to  30  kcal/mol.  The  dielectric  boundary,  which 
should  roughly  correspond  to  the  average  vdW  volume  of  the 
bounded  explicit  waters,  ranged  from  w0(dielectric)  =  w  +  1.8  A 
to  wl (dielectric)  =  w  +  2.2  A. 

Solute  Restraints 


where  N  is  the  number  of  solute  atoms,  x0  is  the  origin  of  the 
simulation  volume,  and  kCOG  is  set  to  100  kcal/mol/A2.  The  rotational 
restraint  (RR)  on  the  solute  involves  restricting  rotation  around  each 
coordinate  axis.  For  example,  the  z-axis  restraint  has  the  form: 


RR 


1 

-  arctan 


(14) 


where  I  is  the  moment  of  inertia  tensor  and  60  is  the  reference 
angle,  which  is  equal  to  zero  when  the  solute  is  initially  trans¬ 
formed  to  standard  orientation.  For  the  layer  simulations  in  this 
work,  kRR  is  set  to  100  kcal/mol/radians2. 


Periodic  Boundary  Simulation  Protocol 

The  periodic  boundary  condition  simulations  in  the  work  were 
performed  with  the  PME  method43  as  implemented  in  CHARMM. 
Local  electrostatic  and  vdW  interactions  were  truncated  with  a 
switching  function  from  7  to  8  A.  In  the  protein  simulation,  a  cubic 
volume  with  at  least  10  A  of  space  between  the  initial  conforma¬ 
tion  of  the  solute  and  the  wall  was  employed  leading  to  dimensions 
of  58  X  58  X  58  A,  and  a  total  of  6231  water  molecules.  The  3D 
fast  Fourier  transform  grid  had  dimensions  of  64  X  64  X  64. 
Langevin  dynamics  simulations  with  a  1  ps-1  friction  constant 
were  run  with  an  extended  system  constant  pressure  algorithm 
using  a  reference  pressure  of  1  atm,  a  piston  mass  of  100  amu  and 
a  collision  frequency  of  10  ps-1. 


Results  and  Discussion 


Because  the  simulation  volume  in  our  method  is  fixed  in  space,  it  is 
necessary  to  restrain  solute  translations  in  both  the  sphere  and  layer 
methods  and  solute  rotations  in  the  layer  method.  The  translational 
restraint  on  the  center  of  geometry  (COG)  of  the  solute  has  the  form: 


In  Tables  1  and  2,  the  accuracy  of  the  multigrid  (MG)  approxi¬ 
mation  is  compared  to  the  standard  switching  cutoff  (STD)  ap¬ 
proach  for  two  model  systems:  a  sodium  ion  embedded  in  a  17-A 
sphere  of  TIP3  water  molecules,  and  two  proteins  surrounded  by  a 
10- A  layer  of  water  molecules.  Three  accuracy  measures  were 
used  to  evaluate  the  multigrid  and  cutoff  approximations  to  the 
combined  electrostatic  and  GB  terms:29  energy  error  (EE), 


A  Hybrid  Explicit/Implicit  Solvent  Method  for  Biomolecular  Simulations 


1973 


Table  2.  Accuracy  of  Multigrid  (MG)  and  Standard  Cutoff  (STD)  Approaches  for  Two  Proteins  Embedded 
in  a  10- A  Solvent  Layer:  (a)  Turkey  Ovomucoid  Receptor  (PDB  Identifier:  lOMT),54  (b)  Trypsin 
(PDB  Identifier:  1TNJ).55 


Method 

Cutoff  start  (A) 

Cutoff  stop  (A) 

%  EE 

%  AFE 

%  MFE 

Time  (s) 

(a) 

MG 

8 

8 

0.55 

4.5 

19 

36 

STD 

26 

28 

0.28 

3.8 

23 

224 

(b) 

MG 

8 

8 

0.25 

5.2 

27.3 

88 

STD 

26 

28 

0.095 

8.0 

55.8 

993 

Definitions  of  error  measures  are  described  in  the  text.  Energies  and  forces  consist  of  the  electrostatic  and  GB  terms 
only.  Last  column  indicates  CPU  time  necessary  to  run  100  steps  of  geometry  optimization. 


I E  ~  Eex act  I 

%EE=100%*  |  |  ,  (15) 

I  Inexact  | 

average  force  error  (AFE), 

2/  mr1/2|Ffact  -  F\ 

%AFE  100% 

I*?  xact  p,  (16) 

and  maximum  force  error  (MFE), 

max[mr1/2|FTxact  -  Ft |] 

%MFE  =  100%  *  ,  (17) 

where  m-  is  the  mass  of  atom  i,  E  is  the  electrostatic  +  GB  energy,  Ft 
is  the  total  electrostatic  +  GB  force  acting  on  atom  /,  the  superscript 
“exact”  corresponds  to  the  infinite  cutoff  result,  and  N  is  the  total 
number  of  atoms.  As  can  been  seen  in  Tables  1  and  2,  a  cutoff  larger 
than  20  A  is  required  to  match  the  results  of  the  multigrid  approach 
that  uses  a  local  cutoff  of  8  A.  Table  1  also  suggests  that  force  errors 
in  the  standard  approach  are  still  significant  even  when  the  cutoff 
exceeds  the  simulation  radius  such  that  all  sodium-water  interactions 
are  included.  In  comparison,  the  multigrid  approach  does  a  reasonable 
job  in  adding  back  the  water  dipole- dipole  interactions  even  though 
the  highest  resolution  multigrid  cell  is  larger  than  a  water  molecule. 
Furthermore,  in  Table  2,  a  comparison  of  the  CPU  times  require  to 
perform  100  steps  of  minimization  for  two  methods  that  have  similar 
accuracy  (MG  8  and  STD  26/28)  suggests  that  the  multigrid  method 
is  approximately  10  times  faster  than  the  standard  cutoff  approach.  All 
CPU  times  stated  in  this  work  were  obtained  from  simulations  mn  on 
an  AMD  Athlon™  MP  2200+  computer  using  the  Portland  Group 
Fortran  compiler  (version  4.0).  Finally,  in  Figure  3,  note  that  the 
multigrid  technique  has  linear  scaling  behavior  with  respect  to  system 
size  for  a  series  of  spherically  solvated  Na+  systems  of  increasing 
size.  It  is  interesting  that  the  standard  procedure  with  large  cutoffs 
actually  does  not  show  linear  scaling  behavior  until  —7000  atoms. 
Another  issue  is  how  accuracy  is  affected  by  increasing  the  system 
size.  As  can  be  seen  in  Figure  3,  the  error  for  multigrid  is  relatively  flat 
vs.  system  size.  The  fact  that  it  does  not  perform  better  with  small 
system  sizes  is  a  probably  a  result  of  the  short  exact  evaluation  cutoff 
(8  A)  and  the  grid  artifacts  mentioned  above. 


Next,  we  evaluate  some  of  the  structural  properties  of  a  spher¬ 
ically  constrained  bulk  water  system  using  the  hybrid  method.  In 
Figure  4,  one  can  see  that  the  hybrid  method  provides  a  reasonably 
homogeneous  radial  distribution  of  oxygen  centers  for  three  dif- 


a) 


Figure  3.  Comparison  of  multigrid  ( a  =  8  A)  vs.  standard  cutoff 
(26/28  A)  for  spherical  solvated  Na+  systems  of  various  sizes,  (a) 
Average  force  error  and  (b)  CPU  time.  Solid  line:  multigrid,  dashed 
line:  26/28  A  switched  cutoff. 
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Figure  4.  Radial  density  profiles  for  spherical  clusters  of  pure  solvent. 
Solid  line:  w  =  12  A,  dashed  line:  w  =  10  A,  dotted  line:  w  =  8  A. 


ferent  system  sizes.  However,  two  notable  deviations  from  homo¬ 
geneity  exist.  First,  oscillations  with  a  magnitude  of  —0.05  and 
periodicity  of  —2-3  A  exist  in  all  three  curves.10  Likely,  this  is  due 
to  edge  effects  that  propagate  from  the  boundary.  The  other  devi¬ 
ation  from  homogeneity  is  the  systematic  decline  of  density  from 
the  center  to  the  boundary.  This  effect  is  probably  the  result  of  a 
moderate  imbalance  between  the  intramolecular  and  boundary 
forces.  A  slight  reduction  in  the  density  of  the  system  reduces  the 
oscillatory  behavior;  however,  the  systematic  density  decline  be¬ 
comes  more  pronounced. 

Figure  5  illustrates  that  the  angular  distribution  of  the  water 
molecules  near  the  surface  deviate  considerably  from  ideal  bulk 
isotropic  behavior,  which  is  a  sine  function,  (1/2)  sin(0).  Most 
importantly,  the  dipole  distribution  graph  is  significantly  skewed 


Figure  5.  Angular  distributions  of  water  molecules  that  are  1  A  or  less 
from  the  surface  boundary  for  a  spherical  cluster  of  pure  solvent  (w  = 
12  A).  Solid  line:  distribution  of  dipole  vectors,  dashed  line:  distribu¬ 
tion  of  normals  to  the  plane  of  water  molecule,  dotted  line:  ideal  curve 
for  both  distributions,  1/2  sin  (9). 


Figure  6.  Radial  distribution  functions  of  oxygen  atom  for  (a)  sodium 
and  (b)  chloride  ions  in  water.  Lines  labeled  “w  =  .  .  .”  correspond  to 
sperical  cluster  calculations. 


towards  75  degrees,  which  roughly  corresponds  to  surface  water 
molecules  with  one  hydrogen  atom  pointing  outwards  from  the 
surface.  The  distribution  of  the  normal  perpendicular  to  the  plane 
is  somewhat  better,  although  significant  deviation  from  isotropy 
exists.  We  tested  angular  correction  terms  for  spherical  boundaries 
similar  to  those  found  in  other  works10,18  and  found  improved 
isotropy  for  both  angular  distributions  (results  not  shown).  How¬ 
ever,  we  did  not  implement  the  correction  terms  for  this  work,  for 
two  reasons.  First,  we  found  that  despite  improving  the  angular 
distributions,  the  dipolar  artifact  near  the  boundary,  which  can  lead 
to  spurious  forces  on  the  solute,  did  not  significantly  diminish.44 
Second,  we  were  unsure  that  applying  an  empirical  angular  cor¬ 
rection  derived  from  spherical  calculations  would  generally  be 
applicable  to  a  layer  where  variations  in  curvature  exist.  Perhaps 
the  angular  distribution  protocol  in  the  SCAAS  method2  could  be 
incorporated,  where  the  desired  distribution  could  be  imposed 
through  angular  restraints. 

Figure  6  shows  the  calculated  radial  distributions  of  the  water 
oxygen  atoms  for  two  simple  monatomic  solute  systems,  Na+  and 
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Table  3.  Charging  Free  Energies  (kcal/mol)  of  Atomic  Ions,  Na+  and 
Cl  ,  in  Spherical  Clusters  of  Different  Sizes. 


Radius  (A) 

Na+a 

Na+b 

cr a 

cr b 

8 

-105.0 

-105.2 

-82.6 

-82.4 

10 

-105.7 

-106.0 

-82.9 

-84.2 

12 

-106.0 

-106.2 

-83.7 

-83.8 

Charging  free  energies  were  determined  by  thermodynamic  integration56 
over  10  equally  spaced  windows  with  charging  parameters,  A  =  0.05, 
0.15,  ...  ,  0.95.  Each  window  consisted  of  20  ps  of  equilibration  and  80  ps 
of  production.  Integration  was  performed  via  a  simple  linear  fit  of  the 
energy  derivatives.  Convergence  errors  are  estimated  to  be  —0.5  kcal/mol 
based  on  a  few  repeated  simulations  that  were  started  with  different 
random  seeds. 

aMultigrid  algorithm  was  used  with  local  cutoff  a  =  8  A. 
bNo  cutoffs  were  used. 


Cl-  ions.  Simulations  of  spheres  with  radii  of  8,  10,  and  12  A  were 
compared  with  PME  calculations  using  a  box  size  of  18  A.  The 
spherical  clusters  are  in  good  agreement  with  PME  with  respect  to 
the  locations  of  the  first  two  peaks,  although  the  cluster  methods 
tend  to  overestimate  the  size  of  the  peaks.  This  problem  is  likely 
related  to  the  results  in  Figure  4,  which  suggest  that  the  water 
density  is  slightly  overestimated  near  the  center  of  the  sphere. 

Table  3  shows  the  calculation  of  charging  free  energies  for 
single  sodium  and  chloride  ions  in  spherical  clusters  of  water 
molecules.  As  can  be  seen  in  this  table,  the  charging  free  energies 
are  approximately  converged  at  around  10  A.41  Furthermore,  ex¬ 
cept  for  one  outlier,  the  multigrid  algorithm  provides  results  close 
to  the  no-cutoff  limit.  In  comparison,  Lee  and  Warshel8  have 
shown  that  simple  radial  truncation  of  electrostatics  without  an 
approximate  long-range  treatment  can  lead  to  charging  free  energy 
errors  greater  than  10  kcal/mol. 

Next,  we  performed  2-ns  molecular  dynamics  simulations  of 
protein  L  in  four  different  environments:  a  layer  without  the  GB 
reaction  field,  a  hybrid  layer,  a  hybrid  sphere,  and  a  PME  cube. 
Figure  7  shows  that  the  root-mean- square  deviation  (RMSD)  tra¬ 
jectories  from  the  original  X-ray  structure  are  similar  between  the 
different  types  of  simulations.  Both  of  the  layer  calculations  seem 
to  wander  less  from  the  X-ray  structure  than  the  sphere  and  PME 
simulations.  The  layer  w/o  GB  simulation  is  especially  restricted 
throughout  most  of  the  2  ns.  The  RMSD  trajectory  results  are 
consistent  with  the  5-factor  data  presented  in  Figure  8  and  Table 
4.  As  can  be  seen  in  Table  4,  the  average  5-factor  calculated  from 
the  layer  w/o  GB  simulation  is  noticeably  smaller  than  the  PME 
calculation,  while  the  sphere  simulation  is  very  close.  The  reduced 
fluctuations  of  the  layer  simulations  may  be  due  in  part  to  the 
angular  restraint  imposed  on  the  moment  of  inertia  of  the  solute 
[eq.  (14)].  Another  possibility  is  that  water  molecules  in  the  layer 
may  have  restricted  flow  leading  to  slight  constraints  on  the 
flexibility  of  the  solute.  Perhaps,  the  layer  w/o  GB  simulation  is 
restricted  more  than  the  hybrid  layer  calculation  because  charge- 
charge  interactions  within  the  solute  are  not  sufficiently  shielded. 
Although  there  is  good  agreement  among  the  four  protocols  for 
calculating  5-factors,  none  of  the  simulation  methods  produces 
5-values  in  quantitative  agreement  with  experimental  data.  This 


discrepancy  can  be  attributed  to  a  variety  of  issues,  such  as  a  lack 
of  anisotropic  components  in  the  fluctuation  calculations,  crystal 
packing  effects,  and  forcefield  errors. 

In  Figure  9,  the  distances  between  the  a-carbons  of  the  final 
simulation  structures  and  the  X-ray  structure  are  presented.  One 
can  see  that  the  hybrid  layer  simulation  has  larger  deviations  than 
the  other  simulation  methods  in  the  main  helical  region  (residues 
23-34)  of  the  polypeptide.  Upon  inspection  of  the  structure,  the 
helix  in  the  layer  simulation  is  still  intact  but  is  shifted  slightly 
relative  to  the  rest  of  the  protein.  One  possible  reason  for  this  result 
is  that  the  helix  dipole  is  interacting  with  the  artificial  dipole  that 
exists  at  the  boundary  of  the  simulation  volume. 

Finally,  the  wall  clock  times  for  the  four  different  protein 
simulations  are  presented  in  Table  5.  As  one  can  see,  the  layer 
protocol  is  almost  three  times  faster  than  the  PME  method  for  this 
protein.  However,  this  is  due  in  large  part  to  the  fact  that  the 
number  of  water  molecules  has  been  reduced  by  a  factor  of  4.5. 
The  computational  times  of  the  layer  and  sphere  methods  suggest 
worse  than  linear  scaling  versus  the  number  of  water  molecules. 
This  result  is  likely  due  to  the  fact  that  the  shape  of  the  simulation 
volume  affects  the  prefactor  of  scaling.  One  reason  why  the  hybrid 
algorithm  is  slower  than  PME  on  a  per  atom  basis  is  because  the 
hybrid  method  must  compute  two  pairwise  energy  terms:  electro¬ 
static  and  GB. 

Our  single  snapshot  error  analysis  results  are  in  good  agree¬ 
ment  with  the  work  of  Skeel  et  al.29  after  one  accounts  for  the  fact 
that,  in  contrast  to  their  protocol,  we  have  used  standard  intramo¬ 
lecular  nonbonded  exclusion  lists.  These  exclusions  remove,  for 
example,  the  large  intramolecular  O-H  electrostatic  interactions  of 
the  TIP3  water  molecule,  which  can  make  percentage  errors  appear 
smaller  by  roughly  an  order  of  magnitude.  The  results  of  Skeel  et 
al.29  suggest  that  the  multigrid  method  does  not  achieve  as  high  a 
degree  of  accuracy  as  fast  multipole  methods.  The  accuracy  in  the 
multigrid  approach  is  mainly  improved  through  extension  of  the 
short-range  cutoff,29  which  can  quickly  become  computationally 
inefficient  because  the  short-range  atom  pair  list  increases  as  the 
cube  of  the  cutoff  radius.  Fast  multipole  methods  have  been  shown 
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Figure  7.  Trajectories  of  the  Ca  RMSD  from  the  X-ray  structure  for 
the  three  types  of  protein  L  simulations,  (a)  Layer,  (b)  layer  (w/o  GB), 
(c)  sphere,  (d)  PME. 
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Figure  8.  Calculated  5-factors  (A2)  of  each  residue  for  the  three  types 
of  protein  L  simulations  compared  to  the  experimental  values.  Calcu¬ 
lated  values  were  obtained  from  the  second  half  of  the  2-ns 
simulations. 


to  provide  highly  accurate  long-range  treatments,  although  some 
difficulties  exist  in  the  evaluation  of  forces  unless  high  precision 
tolerances  are  employed.  The  multigrid  method  has  an  advantage 
over  these  multipole  methods  because  it  utilizes  a  smoothly  vary¬ 
ing  potential  that  leads  to  relatively  stable  molecular  dynamics 
trajectories.  This  smoothness  feature  is  also  useful  if  multiple  time 
stepping  is  desired.29 

An  alternative  approach  to  the  treatment  of  the  long-range 
electrostatic  interactions  is  the  local  reaction  field  (LRF) 
method  of  Warshel  and  co workers. 8  The  LRF  method  builds  a 
very  accurate  long-range  treatment  at  specified  time  intervals  in 
a  molecular  dynamics  calculation.  It  then  uses  the  same  long- 
range  potential  to  rapidly  update  the  long-range  energy  and 
forces  within  each  interval.  It  is  not  obvious  how  the  energies 
and  forces  deviate  compared  to  the  exact  limit  as  the  system 
evolves  from  the  initial  long-range  potential  update  step.  How¬ 
ever,  free  energy  studies  using  the  LRF  technique  have  pro¬ 
vided  accurate  results  vs.  exact  treatments.7,8  Unlike  the  LRF 
method,  the  full  multigrid  quantity  in  our  procedure  is  updated 
at  every  time  step,  such  that  there  are  no  assumptions  regarding 
the  smoothness  of  the  long-range  potential  over  time.8  Instead, 
the  multigrid  approach  induces  a  smoothness  assumption  in  the 
spatial  dimensions  that  likely  impairs  its  ability  to  achieve  high 


Table  4.  Average  Per-Residue  B -Factors  for  the  Last  1  ns  of  Each  of  the 
Four  Different  Types  of  MD  Simulations  of  Protein  L. 


Method 

Average  B -factor  (A2) 

Layer  (w/o  GB) 

16.1 

Layer 

18.2 

Sphere 

19.3 

PME 

19.7 

X-ray 

24.7 

Figure  9.  Ca  distances  to  the  X-ray  stmcture  for  the  final  structures 
generated  with  the  four  types  of  protein  L  simulations. 


accuracy.  In  any  case,  the  most  crucial  aspect  of  the  multigrid 
approach  for  the  purpose  of  this  work  is  that  it  is  the  only 
current  method  that  can  easily  be  adapted  to  handle  the  complex 
GB  potential  term. 

Although  the  initial  results  presented  do  not  comprise  a  com¬ 
prehensive  analysis  of  the  hybrid  method  compared  to  established 
methods,  such  as  PME,  we  would  like  to  distinguish  the  types  of 
applications  where  the  method  may  be  appropriate  and  computa¬ 
tionally  efficient.  First,  the  hybrid  approach  would  be  useful  for 
simulating  very  large  solutes  where  the  box  size  required  to 
encapsulate  the  solute  in  water  molecules  is  prohibitive.11  An 
important  caveat  is  that  the  layer  method  may  not  be  reliable  for 
protein  simulations  where  major  conformation  rearrangements  are 
expected  to  occur,  as  the  simulation  volume  is  held  fixed.  For 
simulations  of  large  flexible  movements  of  molecules  to  be  feasi¬ 
ble,  one  would  need  to  make  intelligent  additions  of  volume  at  the 
beginning  of  the  simulation  or  multiple  restarts  of  the  simulation  as 
the  solute  conformation  changes.  Furthermore,  the  fixed  layer 
approach  would  be  counter-productive  for  folding/unfolding  stud¬ 
ies.45  The  layer  method  may  also  exaggerate  protein  stability 
somewhat  as  was  seen  in  the  results  presented  in  this  work.  In 
contrast,  the  spherical  cluster  approach  should  work  fine  for  cases 
where  large  rearrangements  are  expected  to  occur.  The  sphere 
method  is  about  twice  as  slow  as  the  layer  approach  for  the  protein 
studied  in  this  work.  This  may  be  an  acceptable  trade-off  given  that 


Table  5.  Wall  Clock  Times  for  the  Four  Different  Types  of  MD 
Simulations  of  Protein  L. 


Method 

Time  (days) 

Number  of  water  molecules 

Layer  (w/o  GB) 

3.4 

1378 

Layer 

7.2 

1378 

Sphere 

13.2 

3091 

PME 

19.6 

6231 

For  simulation  details,  see  Methods  section. 


For  simulation  details,  see  the  Methods  section. 
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the  spherical  approach  has  somewhat  better  agreement  with  PME. 
However,  for  certain  biological  systems  that  are  highly  nonspheri- 
cal  (e.g.,  a  DNA  strand),  the  layer  approach  would  likely  be  more 
appropriate  since  the  computational  savings  would  be  more  sub¬ 
stantial.10 

Finite  cluster  and  non-Ewald  PBC  approaches  are  often 
preferred  for  certain  applications  where  complications  regard¬ 
ing  real-virtual  interactions  of  the  solute  may  be  an  issue.46,7 
For  example,  in  the  calculation  of  charging  free  energies  of 
polyatomic  solutes  in  solution,  certain  corrections  to  the  Ewald 
formulation  must  be  employed  to  remove  real-virtual  interac¬ 
tions.47  Other  applications  where  complications  might  arise 
include  the  calculation  of  mutational  free  energies  where  net 
charge  changes  occur  and  the  computation  of  pKas  of  ionizable 
residues.  Many  large-scale  charging  calculations  to  date  have 
been  performed  with  non-Ewald  PBC.48,49  Unfortunately,  non- 
Ewald  PBC  methods  require  relatively  large  cutoffs,  and  hence, 
increased  computational  effort  for  the  nonbonded  electrostatic 
terms  vs.  PME.  The  hybrid  method,  on  the  other  hand,  thanks  to 
the  multigrid  enhancement  and  reduction  of  explicit  water  mol¬ 
ecules,  achieves  similar  computational  efficiencies  as  PME 
without  incurring  the  artificial  periodicity  of  the  solute.  In  fact, 
in  a  separate  article,41  we  applied  the  hybrid  approach  towards 
calculating  solvation  free  energies  of  fixed  protein  conforma¬ 
tions.  The  resultant  data  are  then  used  to  benchmark  the  accu¬ 
racy  of  various  Poisson  implicit  solvent  models. 

Similar  to  other  finite-cluster  methods,  such  as  GSBP  and 
SCAAS,  the  hybrid  method  is  not  limited  to  explicit  solvation 
of  the  entire  solute.20,50  This  is  a  major  computational  advan¬ 
tage  vs.  PBC  approaches.  For  example,  in  certain  applications, 
such  as  protein-ligand  binding,  protein-protein  binding,  and 
mutational  studies,  the  simulation  volume  can  be  defined  in 
terms  of  a  subset  of  “active”  atoms.  The  nonactive  atoms  can 
then  be  treated  with  a  regular  GB  model.  Because  the  simula¬ 
tion  volume  is  fixed  in  the  current  procedure,  the  nonactive 
atoms  should  be  fixed  also.20 

In  principle,  the  hybrid  method  can  be  extended  to  allow  for 
a  dynamically  adjusting  layer  that  can  achieve  constant  pres¬ 
sure.10,11  There  are,  however,  several  difficulties  with  a  con¬ 
stant  pressure  approach  that  may  need  to  be  addressed.  First, 
large  conformational  changes  of  the  solute  may  lead  to  signif¬ 
icant  changes  in  the  geometric  volume  of  the  simulation  vol¬ 
ume.  Thus,  it  may  be  necessary  to  add  or  delete  water  molecules 
using  a  constant  chemical  potential  approach  such  as  a  grand 
canonical  Monte  Carlo  procedure.51  Moreover,  a  dynamic  layer 
adjusting  to  local  pressures  may  squeeze  out  waters  from  un¬ 
favorable  locations  near  a  hydrophobic  region  of  the  solute.11 
This  problem  suggests  that  the  layer  widths  for  each  heavy  atom 
may  have  to  be  restrained  from  becoming  too  small.  In  addi¬ 
tion,  coupling  of  the  boundary  to  the  solute  positions  may 
lead  to  decreased  fluctuations  of  the  solute  atoms  due  to  the  fact 
that  waters  near  the  boundary  will  be  dragged  around  as  the 
solute  moves.  Finally,  the  simplicity  and  computational  effi¬ 
ciency  of  using  a  Born  radii  grid  to  interpolate  atomic  Born 
radii  will  have  to  be  sacrificed  if  the  boundary  needs  to  be 
updated  often. 


Conclusion 

We  have  introduced  a  new  algorithm  for  performing  biomolecular 
simulations  of  a  solute  surrounded  by  an  irregularly  shaped  layer 
of  water  molecules.  Unlike  previous  layer-based  approaches, 
proper  reaction  field  treatment  is  achieved  by  encapsulating  the 
explicit  simulation  volume  in  an  implicit  solvent  described  by  the 
GB  theory.  Significant  computational  enhancement  of  this  ap¬ 
proach  is  achieved  through  the  use  of  a  pairwise  multigrid  tech¬ 
nique  that  has  been  extended  to  incorporate  the  GB  model.  This 
method  overall  can  be  faster  than  the  conventional  PME  technique 
in  cases  where  the  number  of  water  molecules  can  be  significantly 
reduced.  We  foresee  the  application  of  the  hybrid  technique  in 
simulations  where  a  finite  cluster  of  water  molecules  is  preferred, 
for  example,  the  calculation  of  binding  affinities  of  protein-ligand 
and  protein-protein  complexes,52  and  the  computation  of  muta¬ 
tional  free  energies.53 
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